clear; clc; close all;

k = 3;
N = 1000; % 网格分辨率

mu_list = linspace(-pi, pi, N);
eta_list = linspace(-pi, pi, N);
[MU, ETA] = meshgrid(mu_list, eta_list);

TYPE_STR = strings(N, N);

parfor i = 1:N
    temp_STR = strings(1, N); % parfor 里避免直接写二维变量，先临时存储
    for j = 1:N
        mu = MU(i,j);
        eta = ETA(i,j);
        a1 = 1;
        a2 = 2- k * cos(eta)- k * cos(mu);
        a3 = (2 - k * cos(mu)) * (1 - k * cos(eta));
        rts = roots([1 a1 a2 a3]);
        temp_STR(j) = Classify(rts);
    end
    TYPE_STR(i, :) = temp_STR; % 一次性赋值，避免数据冲突
end

% 映射字符串到数字（方便绘图）
labels = ["USP (2P1N)", "USP (1P2N)", "USF (1P2N)", "USF (2P1N)", "SNF (0P3N)"];
TYPE_NUM = zeros(N, N);
for i = 1:N
    for j = 1:N
        idx = find(labels == TYPE_STR(i,j));
        if isempty(idx)
            TYPE_NUM(i,j) = 1; % Other映射为1
        else
            TYPE_NUM(i,j) = idx;
        end
    end
end

% 自定义颜色表
cmap = [
    1 0 0;        % USP (2P1N) 红
    0 0 1;        % USP (1P2N) 蓝
    0 0.7 0;      % USF (1P2N) 绿
    0.5 0 0.5;    % USF (2P1N) 紫
    0.7 0.7 0.7]; % SNF (0P3N) 灰

figure; hold on;
imagesc(mu_list, eta_list, TYPE_NUM);
set(gca, 'YDir', 'normal');
colormap(cmap);
colorbar('Ticks', 1:length(labels), 'TickLabels', labels);

caxis([0.5 5.5]); % 限制颜色条范围

xlabel('\mu');
ylabel('\eta');
title('Fig.1 稳定性分类图');
grid on;

xlim([-pi, pi]);
ylim([-pi, pi]);





